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Summary 

\ 

A slowly convergent stationary iterative process can be accelerated by explicitly annihilating 
(i.e., eliminating) the dominant eigenvector component of the error. The dominant eigenvalue 
or complex pair of eigenvalues can be estimated from the solution during the iteration. The 
corresponding eigenvector or complex pair of eigenvectors can then be annihilated by applying 
an explicit Richardson process over the basic iterative method. This can be done entirely in 
real arithmetic by analytically combining the complex conjugate annihilation steps. We illustrate 
*by applying the technique to an implicit algorithm for the calculation of two-dimensional steady 
transonic flow over a circular cylinder using the equations of compressible inviscid gas dynamics. 
This demonstrates the use of explicit annihilation on a nonlinear problem. 


Introduction 

A great deal of work is being done to develop efficient algorithms for steady-state problems in 
computational fluid dynamics. The most widely used class of algorithms for steady problems consists 
of time-like methods that are marched to a steady state. Our purpose here is to point out a simple 
and easily applied technique from linear algebra which can in some cases produce a remarkable 
speedup in algorithms that march to the steady state. The idea may be briefly summarized as 
follows. In the later stages of the marching algorithm, the iteration may be behaving linearly. If so, 
one can estimate the dominant eigenvalue of the underlying iteration matrix and “annihilate” the 
eigenvector corresponding to this eigenvalue by a special iteration step. Even when the eigenvalue 
and eigenvector are complex, all the computations can be performed in real arithmetic. 

Some notation is developed and some facts from linear algebra are reviewed in the next section. 
A particular algorithm for the Euler equations is presented in the third section, and an application 
of annihilation to steady transonic flow over a circular cylinder is described in the final section. We 
wish to emphasize that the application of annihilation is not restricted to the particular marching 
algorithm then we use for illustrative purposes. 

We extend our thanks to Harvard Lomax for introducing us to the idea of eigenvector annihilation. 

Review of Linear Algebra 

In this section we establish our notation and review some elementary facts from numerical linear 
algebra. To begin, consider the matrix problem Ax = b where A is nonsingular. Let A = M — N 
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be a splitting of A, where M is nonsingular, and consider the iteration 


= Nx^ -f 6, n > 0, j^given, 

or M/^x^ = —Ax^ + 6, (1) 

and := x” + Ax^. 


(The notation “a := 6” is used here to mean that a is defined as b.) 

Let x* satisfy Ax* = 6, and denote the error by e'* := x* — x". Then and so the 

iteration converges if and only if the spectral radius of N is less than 1, p(M~^ N) < 1. 

Suppose now that N is diagonalizable with eigenvalues X,- and corresponding eigenvectors v,*. 
Take |Xi| > IX 2 I > .•• > \^m\, and write o,V|. Then we clearly have 


e"+* = 

t==l 

m 

= X,(a,vi + ^ a.(Xf/Xi)*'t',). 

1=2 


( 2 ) 


Evidently the component of the error in the direction of Vi is the slowest decreasing component (if 

|Xi| > (Xal)- 

The component of the error in any direction Vj- can be annihilated as follows. ARer computing 
Ax”"*"*, define x"+*+* := + aAx"'^^, where cr is as yet arbitrary. (We call this a Richardson 

step with parameter a.) Then we have 

gn+*+l gn+fc _ ^^j-n+fc 

= e"+* - + 6) 


= [/-<rM-'(M-fV)]e"+'' 

I 

= “i-XJKI ~(r) + <rX,]vi. 

t 


Thus if we choose a = 1/(1 — Xy), the component of in the direction of vj is zero. 

Any component vj can in theory be annihilated, but in practice the important component is the 
eigenvector v\ associated with the dominant eigenvalue Xi. 

Of course, we do not know Xi, but we can use the sequence Ax^ to estimate it, since Ax'^'^^ = 
M~^ N Ax^ . Estimating Xi from the sequence Ax'^ is simply using the power method to find the 
dominant eigenvalue of a matrix. Thus we estimate Xi by Xi (Ax”'*"^)r/(Ax^'*"^”^ )r for some 
appropriate component of the update Ax. (The residual := —Ax^~^^ + 6 could also be used 
for the estimation, since = Ae^'^^ .) In summary, we plan to estimate Xi from the sequence 
Ax^ and take a Richardson step with parameter cr = 1/(1 — Xi ). 

The general condition for stability of the Richardson step with parameter a is |l — (j( 1 — X,)| < 
1 for all eigenvalues X,- of N . This ensures that no component of the error is amplified 
by the Richardson step. It may well happen that the Richardson step we propose is “unstable'’ 
in the senvSe that some components of the error are magnified. We note, however, that in the 
later stages of an iteration the subdominant components of the error are liable to be so small 
that a certain amplification of them is acceptable when accompanied by a large decreaj>e in the 
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dominant component of the error. Also, further ordinary steps will rapidly reduce these subdominant 
components again. To be more specific, suppose = aiVi + ev,-, t ^ 1. Then we easily find 

= Q!i(l — 6T-h crXi)i;i + e(l — (T -f cr\i)vi 

= if <T= 1/(1-Xi). 

1 — Ai 

Thus, the v,- error component is multiplied by a factor of (X,- ~ ^i)/(l — Xi). We expect Xi to be 
close to 1, so this factor could be large, but the factor of e will compensate if it is small enough. 

We remark here that in our application, A is not constant, but A = A(x”). Nevertheless, the 
assumption of linear behavior will be approximately satisfied in the later stages of the iteration. 
We will not attempt annihilation unless the iteration seems to be behaving linearly. An empirical 
criterion for determining this is described in the last section. 

We will now show how annihilation can proceed in real arithmetic even when the dominant 
eigenvalue is complex. It is quite possible that the dominant eigenvalue could be complex, and thus 
{since we will assume that A is real) that there is a dominant pair of complex conjugate eigenvalues, 
say X 2 = Xi, and |Xi | = IX 2 I > |Xs| > . . . > |Xy„|. Then if we write fiiViy we see 

^^n+k-p _ xJ-Vi Vi + X^"^2V2 + OdXsl''-") 


for p = 0, 1,2. Thus, for any real c and d, 

Ai"+* + cAx'*+*-' + rfAz”+*-2 = (X? + cX, + d)vr 

+ X 2 ^ 02{^2 + ^)V2 + *)• 

The coefficients of Vi and V 2 will vanish if c = —2Re\\ and d = |Xi|^, for then X^ + cX + d = 
(X-X 1 KX-X 2 ). 

In practice, we can pick indices i j and find c and d such that 

(Ai"+*),- + c(Ax"+*-* ), + d(Ax"+''-^)( = 0 
(Ax"+*)y + c{Ax"+*-* )y + d(Ax"+''-2)y = 0. 

This is a linear system for c and d. Once c and d are known, Xi = —c/2 + iy/ d — c^/A and two 
Richardson steps can be performed by defining a := 1/(1 — Xi) and putting 

j.n+fc+1 j.n+fc ^ 

j.n+Ar+2 ^n+fc+1 ^ ^ ( ) 


where Ax is still given by one step of the original iterative scheme, namely, Ax”"'"*' := 

(— /Ix”"^* + 6) and Ax"’*'*"*"' := (— + 6). By the calculations done previously, 

this will eliminate the Vi and V 2 components of the error. Furthermore, there is no need for complex 
arithmetic, for the computation can be reorganized as follows: 

A/(x"+''+2 - x’*+*) = A/(9 Ax"+''+* + <tAx'*+*) 

= a(-/lx"+''+' + 6) + <7(-ylx"+* + 6) 

= a[-A(x”+'‘ + (tAx"+'')]? 6 - aAx"+'' + ab 
= —{a + a)/4x"‘*‘* + (<r + c)b 

-1<t|MM-*(-/1x"+* + 6) 

= 2Re<7 r"+'' - \a\^ AM~^ . 


I 
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Thus, two complex Richardson steps can be performed entirely in real arithmetic, at a cost of 
(1) solving two linear systems with coefficient matrix M and (2) performing one matrix-vector 

multiplication with coefficient matrix A, This is the same as the cost of two steps of the basic 
iteration process. As for stability, a component of in the direction of the vector vj, j > 2, is 
multiplied by the factor _ 

(x.-xi)(xy-y.) 

(I-Xi)(l-Xi) • 

Again, this factor could be large, but we expect to apply the annihilation in the later stage of the 
iteration when the coefficients of vj are very small and some amplification would be tolerable. 

We might remark that the presence of a dominant complex conjugate pair of eigenvalues is signaled 
by oscillations in the norm of the residual. For example, if = aiAv\ H- a 2 Av 2 = 2Re(aiAvi), 
then writing Xi = we see that in any norm 

Thus, the vector aiAvi is rotated through an angle k9 and projected on the real axis, which will 
give oscillations. Furthermore, if ^ = 'npjq in lowest terms, then the period of the oscillation will 
be q. Examples of oscillatory decreasing residuals (which signal the presence of a dominant complex 
conjugate eigenvalue pair) can be found in references 1-3. 

We close this section with some remarks on the relation of annihilation to previous work. It is easy 
to show that in the case of a scalar sequence (x®, our “annihilation” step is, exactly, an 

Aitken ^^-extrapolation step. Wilkinson (ref. 4) discusses the application of the Aitken 5^ technique 
to the power method for estimating eigenvalues and eigenvectors, where x” is a vector. The formulas 
there are different from the ones we have used: in Wilkinson’s formulation the Aitken acceleration 
step is applied to each component of the vector separately (in effect, an “eigenvalue” is assumed for 
each separate component of the sequence of vectors), whereas in our formulation a single eigenvalue 
is estimated and used for all components of the vector. We are unable to say which strategy is 
better. Furthermore, in Wilkinson’s formulation it is unclear how to handle the case of a dominant 
complex conjugate pair of eigenvalues (it is shown how to estimate the eigenvalues as roots of a 
quadratic equation, but the acceleration of the convergence of the vectors is not discussed). The 
idea of annihilation is elementary but appears to have been neglected. Lyusternik (ref. 5) considered 
a procedure that amounts to annihilation in the case in which the eigenvalues are all real, but he did 
not discuss the complex case. Hyman and Manteuffel (ref. 6) described an algorithm very similar 
to annihilation for accelerating slowly convergent iterative methods. In their algorithm, eigenvalue 
estimates are obtained via a Krylov sequence technique. 

Many authors have been concerned with the harder problem of optimizing relaxation schemes 
(getting all the eigenvalues of the iteration matrix as far inside the unit circle as possible). An 
adaptive procedure for minimizing the norm of the iteration matrix was presented by Manteuffel 
(ref. 7); he estimated complex eigenvalues by the same technique we use above. To reiterate our 
point of view: we do not ask for the “optimum” method in some class of methods, nor do we insist 
that our annihilation steps be stable. We simply estimate eigenvalues and perform annihilation. 
Moreover, we are willing to accept amplification of some error components in the annihilation step; 
these components will be rapidly reduced anyway in subsequent normal iteration steps. 
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An Iterative Method for the Euler Equations 

In this section we will describe an iterative method for some steady-state problems in fluid 
dynamics. The governing equations are the steady two-dimensional Euler equations of compressible 
inviscid gas dynamics, which express the conservation of mass, momentum, and energy. They can 
be written in the form 

where Q is the four-vector Q = (p, pu, pv, e)^. Here p is density, u and v are Cartesian velocity 
components, and e is total energy per unit volume. The functions E,F : D CZ are 

nonlinear functions given by 


^(Q) = {p^< + P> p^^’ “(e + p)) 

F(Q) = (^pv,puv,pv^ +p,v{e + p)^ . 

The pressure p is defined by p = [7 — l)(e— ^p(u^ -f v^)] where 7 is a constant. Appropriate 
boundary conditions must be adjoined to the differential equation to complete the specification of 
the problem. It is not important for our purposes what these boundary conditions are, so we will 
not discuss them further. In order to handle a curved geometry, we map from the physical (x,p) 
domain to a “computational” (f, 7 ) domain; the computational domain is usually taken to be a 
rectangle, for ease in differencing. In the (^, 7 ) coordinates, the transformed equations retain the 
strong conservation law form of equation ( 11 ) (ref. 8 ), becoming 


dE(Q) dF(Q) 
^ dv 


( 12 ) 


Iterative procedures for solving the steady-state equations (12) often take as their starting point 
the unsteady equations 


dQ dE(Q) dF(Q) 
dt ^ ^ dff 


(13) 


If we choose Euler implicit differencing in time (hoping for good stability properties) we get the 
iteration 




Q =Q 




(14) 


where 6 ^, 6n are spatial difference operators, h := Af, and n denotes a time level. To avoid iteration 
on the nonlinear flux terms, E^^^ and F on the right-hand side of equation (H) are expanded 
about values at level n using a Taylor series 


£"*■ _ E* + ^ 
dQ 

A n * n + 1 - w 

-Q 


AQ" + (^((AQ")") 


The locally linearized form of equation (14) can then be written in “delta” form as 


(15) 


(/ + A(5(A" + 5,b")]AQ'* = -h(S(E’' + «,f"), 


(16) 
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where and are the Jacobian matrices dE and dP^jdQ'^ , respectively. 

Each step of equation (16) involves the inversion of a large block-banded matrix, with half- 
bandwidth equal to the number of points in one direction of the mesh. The amount of work required 
for this is unacceptably large; instead the matrLx is approximately factored (ref. 9) For example, if 
and bf! were three-point central difference operators, the left-hand side of (16) might be factored 
via 

[/ + hib^A^ + 5 ^ b ”)1 = (/ + hb^A")(I + hbr^B'') + O(h^) ( 17 ) 


as the product of two block-tridiagonal matrices. 

The calculations we will show come from an algorithm that uses flux-vector splitting (refs. 10, 
11), which is based on separating positive and negative characteristic directions to allow the use of 
one-sided spatial difference operators. The Euler equations (13) are hyperbolic, and E = AQ can 
be decomposed as 


E = AQ = XAX-^Q 
= X(A++A-)X-‘Q 
= {A^ +A~)Q 


(18) 


= E'^ + E , 


where A is the diagonal matrix of eigenvalues of A, and and A” are the positive and negative 
parts of A, respectively. F can be similarly decomposed. Defining A^ := DE^fdQ, and similarly 
for A , B^ j and B , the unfactored scheme is 


\I + h(6\A^ + + 6{,B )”]AQ’‘ 

= -h{d\E'^ + 6^(E~ + 


where b^ and b^ are backward and forward spatial difference operators. On the left-hand side, 
the forward difference operators are separated from the backward, resulting in an approximate 
factorization into lower and upper block-triangular matrices. First-order spatial differences and 
are used on the left-hand side and second-order differences b^ and b^ are used on the right-hand 
side, to give second-order accuracy in the steady state. We have, for instance. 


A* A,- := A,- - A,_i 

6»^Ei := (3Ei - 4E._i + E,_2)/2. 


The full scheme is thus 

[/ + h(A\A'^ + + /i(A^A" + AJ;b")"]AQ" 

= -h(6\E'^ + 6{E~ + 

In the notation discussed in the preceding section. 


A/= + + A‘c'^)'*)(/ + /i(A^.4 +A{B )"] 

A = + 6'’^B'*' + S{,B~)’'. 


It is clear that A is not a constant matrix, but if we wait until the later stages of the iteration, A will 
be changing very slowly, and we should be able to apply annihilation with a fair degree of success. 
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Application to Steady TVansonic Flow 

The annihilation idea has been applied to the numerical method given in the preceding section 
for the Euler equations. In this section we will describe the physical problem, show how to apply 
annihilation, and give some results showing the success of annihilation. We will also mention some 
questions and problem areas that our experiences have identified. 

We performed calculations for fiow over a circular cylinder with symmetry imposed between the 
top and bottom, so that only the top half of the region was calculated. Figure 1 shows the grid 
used, including the two points below the symmetry line in front of and behind the cylinder used to 
impose the symmetry condition. At the far-field boundary, 16 diameters away from the cylinder, 
all flow variables are fixed at their free-stream values. Boundary conditions at the body consist of 
setting the normal velocity to zero, taking surface density and tangential velocity from the point 
above, and calculating pressure from a conservation of momentum relation. 

At a free-stream Mach number of 0.5, a shock forms on the cylinder. A steady solution for this 
case is shown in figure 2. The shock has introduced rotationality into the flow and caused inviscid 
flow separation on the back side of the cylinder. 

The application of eigenvector annihilation to the iterative method (eq. (20)) requires some ad 
hoc decisions on how to estimate the dominant eigenvalue pair. The method presented here is one 
of many strategies that could be envisioned. As described in the second section (Review of Linear 
Algebra), one component of AQ is needed at two points in the field, corresponding to subscripts x 
and j in equation (6), for three consecutive iteration steps. We used the first component (density) 
of AQ. Four grid points were chosen, each a third of the way in from a corner of the computational 
domain (see fig. 1). The two points closest to the cylinder form one pair, the other two another 
pair. In this way we can obtain two estimates of the dominant eigenvalue, from different regions in 
the grid. The two points from each pair are separated in an attempt to minimize any local coupling 
effects. During each iteration^ eigenvalue estimates are made from the two pairs of points. If the 
real and imaginary components of the estimates differ by less than 5%, based on the modulus of the 
first estimate, the two estimates are averaged to produce a candidate X. Annihilation is performed 
if this candidate X is within 5% of the candidate X from the previous iteration. Thus, our criterion 
for linear behavior involves a consistent eigenvalue estimate from widely separated points in the grid 
over four iterations. 

The estimated eigenvalue is in general complex; hence, two Richardson steps are required. An 
outline of the steps performed is given below, using the definitions of the matrices M and A from 
equation (21): 

1. Ag” = M{Q"r^A{Q")Q" 

2. Decide whether to annihilate. If yes, 

3. =Q’* + (l<T2l/(2fte<7)]AQ" 

4. Apply boundary conditions to 

5. = M(g")-M(g'*)g"‘^' 

6. g”^^ = g" + 2/?e<TAg”^* 

- n4-2 

7. Apply boundary conditions to Q 

The result of applying annihilation, using this strategy, to the cylinder problem is shown in figure 
3. In this figure we compare the result of no annihilation with the result of starting the annihilation 
strategy at n = 500. In figure 3 annihilation steps were performed at n = 504, 540, 577, 616, 670, 
704, 869, 892, and 904. The convergence rate over the last 300 iterations for the curve without 
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annihilation was about 0.9958, whereas for the curve with annihilation the convergence rate was 
about 0.9881. We remark that if only one annihilation step is performed, the residual drops sharply 

but eventually resumes converging at a rate of 0.9958. This may be due either to error in the 
estimation of the eigenvalue (so that the dominant eigenvector component of the error was not 
completely annihilated) or to nonlinear feedback effects stemming from the fact that the iteration 
process is truly nonlinear. 

We close by mentioning some questions and problem areas that have arisen during this work. 
First, the choice of the strategy used to estimate the eigenvalues — ours was heuristic — is important. 
Can better ones be devised? In this regard we would like to mention the note by Jones (ref. 12), ho 
gives a statistical criterion based on a serial correlation coefficient for deciding when to employ the 
Aitken technique. The application is to a real scalar sequence; it would be useful to have a similar 
statistical criterion for a trigonometric sequence modulated by a geometrically decaying term, but 
we are unaware of a serial correlation coefficient in this case. 

Second, it is difficult to decide when to annihilation should be started. In figure 4 we show the 
result of starting the annihilation strategy at n = 0. In this case there were 19 annihilation steps: at 
n = 129, 168, 208, 235, 268, 310, 369, 476, 491, 626, 648, 656, 794, 896, 918, 927, 934, 951, and 984. 
Evidently this strategy was successful, but not much more so than that of waiting until n = 500 to 
start annihilating. This may be due to an incorrect estimation of eigenvalues or to the nonlinearity 
of the process. Third, stability questions arise when performing annihilation. We see in figure 4 
some sharp increases in the residual at certain annihilation steps, and annihilation strategies that 
do not allow these jumps might be preferred. Finally, the annihilation procedure requires an extra 
(third) level of computer storage beyond that normally needed for the iterative procedure; this may 
be a problem in cases where extra space is scarce. 
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Fig. 1. Exponentially stretched 42-by-31 grid about a circular cylinder, showing 
points used for eigenvalue estimation. 



Fig. 2. Streamlines and sonic line for steady flow about a circular cylinder at a 
free-stream Mach number of 0.5. 
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Fig. 3. Comparison of convergence histories without annihilation and with annihila- 
tion starting at n = 500. 
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Fig. 4. Convergence history with annihilation starting at n = 0. 
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